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ABSTRACT 

The essentially non-oscillatory (ENO) finite difference scheme is applied to systems 
of conservation laws u t + f(u) I = 0 of mixed hyperbolic-elliptic type. A flux splitting 
f(u) = f + (u) + f (u), with the corresponding Jacobi matrices having real and pos- 

itive/negative eigenvalues, is used. The hyperbolic ENO operator is applied separately on 
f+ ( u )* and on f (u) x . The scheme is numerically tested on the van der Waals equation in 
fluid dynamics. We observe convergence with good resolution to weak solutions for various 
Riemann problems, which are then numerically checked to be admissible as the viscosity- 
capillarity limits. We also observe the interesting phenomena of the shrinking of elliptic 
regions if they are present in the initial conditions. 


1 Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in Science and 
Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. Research also was partially 
supported by NSF Grant DMS-88-10150, NASA Langley Grant NAG-1-1145 and AFOSR Grant 90-0093. 
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1 Introduction 


The system of conservation laws 


U < + f ( U )x = 0 

u(x,0) = u°(x) 


( 1 . 1 ) 


is hyperbolic if the Jacobi matrix has real eigenvalues and a complete set of eigen- 

vectors. In recent years there have been a lot of activity in designing stable and accurate 
numerical methods for solving systems of hyperbolic conservation laws. The ENO (essen- 
tially non-oscillatory) high order finite difference method, [5], [6], [12], [13], is one of the 
successful approaches. The philosophy of ENO schemes is to use upwinding and adaptive 
stencils, based on the local “wind” direction (the sign of the relevant eigenvalue) and the 
local smoothness, in each of the local characteristic fields. ENO schemes can resolve the 
interactions of shocks and waves and complicated wave structures quite well, according to 
the numerical examples for the hyperbolic Euler equations of a polytropic gas in [6], [13]. 
ENO schemes are formally high order accurate, measured by local truncation errors. In some 
cases the actual order of accuracy may degenerate due to the frequent switching of stencils in 
the linearly unstable regions [9], [14]. There is a simple modification of ENO schemes, with 
no additional computational cost, [14], which overcomes this accuracy degeneracy problem. 
However, our experience seems to show that the original ENO scheme has very good resolu- 
tion for nonlinear systems [13], [14], suggesting that the modification may not be necessary 
in such cases. 

If the Jacobi matrix in (1.1) has complex eigenvalues, the system becomes elliptic. 
The initial value problem (1.1) may not be well posed in the elliptic regions. However, in 
applications, (1.1) may still be of mixed hyperbolic-elliptic type. Examples include equations 
in fluid dynamics [17], elasticity [8] and the partial differential equations related to Lorenz 
systems [7], just to name a few. In this paper we use the van der Waals equation in fluid 
dynamics 


vt + p(w)x = 0, w t - v x = 0 
u(x,0) = v°(x), u>(x,0) = tu°(x) 


( 1 . 2 ) 


with 


p(w) = 


RT 


a 


(1.3) 


w — b w 2 

where R,T,a,b are all positive constants, for our numerical examples. See, for example, [17] 
for details. Equation (1.2) corresponds to (1.1) with u = ( v,w ) T , f(u) = (p(w),—v) T . The 
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two eigenvalues of the Jacobi matrix are ±y— For an ideal gas, the pressure p(w) 
is a decreasing function of w, resulting in a hyperbolic system (1.2). However, during the 
co-existence of gas and liquid, p'(w) may become positive within an interval, as in the case 
(1.3) with suitable parameters (Figure 1), making the system (1.2) elliptic in this region. The 
mathematical ill-posedness of the system in the elliptic region reveals the physical fact that 
the state in the elliptic region is not stable, and it typically evolves into “phase transitions” , 
i.e., jumps across the elliptic regions in the weak solution. As in the hyperbolic case, there can 
be more them one weak solution, and there are efforts in the literature, e.g., [15], [16], [8], [2], 
[3], [10], [11], to establish admissibility criteria with the goal to single out one “physically 
relevant” weak solution. A weak solution u = ( v,w ) of (1.2) is called admissible as the 
viscosity-capillarity limits in [15] if it is the bounded a.e. limit of u e = (v e ,w e ), satisfying 

' v\ + p(w e ) x = ev xx - e 2 Aw xxx , 

< W\-V% = 0 (1.4) 

v € (x, 0) = v°(cc), w € (x, 0) = u/°(a;) 

as e -> 0 + . Here A is a positive constant, usually between 0 and 1. The analysis usually 
starts with the Riemann problem 


(„(*,<», »(*,<>)) = 1 x >° 0 


(1.5) 


whose solutions contain most of the rich features (shocks, phase boundaries, rarefaction 
waves, etc.) for more general initial conditions. It is easier to consider each jump (shock, 
phase transition) separately, with a jump called locally admissible if it is the limit of the 
travelling wave solutions of (1.4). A weak solution is then called locally admissible if each 
of its jumps is locally admissible. Admissible or locally admissible solutions of the Riemann 
problems, e.g. [15], [10], [3], typically contain phase boundaries (discontinuities jumping across 
the elliptic regions) but do not achieve values inside the elliptic region. More general initial 
data is more difficult to analyse. It is interesting to know what happens if the initial condition 
contains elliptic regions. i,From a computational point of view, if a shock capturing method, 
such as the ENO method in [6], [13] is to be used, discontinuities are typically spread out in 
two or three points, hence points can sit in elliptic regions even if the exact solution jumps 
across it. 

One of the main ingredients of ENO schemes, and of many other non-oscillatory schemes 
such as TVD (total-variation-diminishing) schemes [4], is the approximation in each of the 
local characteristic fields. If the system becomes elliptic, local characteristic decomposition 
is no longer available. In Section 2 we propose a flux splitting f(u) = f + (u) + f - (u), with 
the corresponding Jacobi matrices having real and positive/negative eigenvalues. This 
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is similar to the flux splitting used for hyperbolic systems, for example the Lax-Friedrichs 
splitting and the van Leer splitting [19]. The hyperbolic ENO operator is then applied 
separately on f + (u) x and on f“(u) x , using the component version in the elliptic regions, 
since the characteristic decompositions of have no physical meanings in such regions. 

This idea is described in the next section. In Section 3 we numerically test the scheme on 
the van der Waals equation (1.2)-(1.3). We observe convergence with good resolution to 
weak solutions for various Riemann problems. We then numerically check that these weak 
solutions are admissible as the viscosity-capillarity limits by computing (1.4) with a sequence 
of decreasing e. We also compute the solution with smooth initial conditions, and observe 
the interesting phenomena of the shrinking of elliptic regions if they are present in the initial 
conditions. 


2 The Numerical Scheme 

We start with a flux splitting 


f(u) = f + (u) + f (u) (2.1) 

with the requirement that the Jacobi matrix has only real and positive eigenvalues, 

and likewise that the Jacobi matrix has only real and negative eigenvalues. If the 

system is hyperbolic, the simplest way to achieve such splitting is due to Lax-Friedrichs 

f ± (u) = ^(f(u) ± aw)) a = I Ai(u) I (2.2) 

where A,(u) are the eigenvalues of the Jacobi matrix . For special classes of hyperbolic 
systems, such as the Euler equations of a polytropic gas, more sophiscated splittings with 
better physical meanings are available, e.g., van Leer’s splitting [19]. For an elliptic system, 
the simple splitting (2.2) no longer works. In fact, any splitting with the Jacobi matrices 
island commuting with each other, as is the case in (2.2), will probably fail, 

because commuting matrices with distinct eigenvalues can be simultaneously diagonalized, 
hence the eigenvalues of their sum are simply the sum of their corresponding eigenvalues. 
However, a splitting similar to (2.2) with the scalar a replaced by a diagonal matrix: 


( 


f± ( U ) = o( f ( U ) ±aU )> a = 


Ot 2 


/ 


(2.3) 
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can usually yield a required splitting. For example, the flux f(u) in the van der Waals 
equation (1.2) can be split successfully using (2.3) with: 


a 2 = max 

w 


max 


o, 


yOkf 2 — 4 p'(iy) — A/' 


OL\ — a 2 H" A/ 


(2.4) 


and 


M = a max y'max^, p'(u>)), o> 2 (2-5) 

The idea is to make an ansatz g(u) = M(v, 0) T , then try to find the smallest possible M 
such that the Jacobi matrices both have real and distinct eigenvalues. This leads 

to M given by (2.5). Once this is done, it is easy to use the Lax-Friedrichs idea, i.e., to add 
and substract au with a suitable a, to accomplish the splitting. 

Equipped with the splitting (2.1), one can then apply any successful hyperbolic approxi- 
mation techniques separately to f + (u) and f“(u). The only exception is that characteristic 
decompositions should not be performed in elliptic regions, since the characteristic directions 
of and af ~^ u ) do not have any physical meaning. In this paper we apply the ENO 

techniques developped in [12], [13] to f + (u) and f~(u). 

We summarize the algorithm briefly in the following. More details can be found in [12], 
[13]. 

(1) The spatial operator 


-f(u) I = -f + (u) x -f (u). 
is approximated by a conservative flux difference 

L(u) j = L+(u) i + L"(u) i 

where j is the grid index (the grids are located at Xj = j Ax), and the resulting ODE 

lH = L(u) > 


(2.6) 


(2.7) 


( 2 . 8 ) 


is discretized in the time variable t by a class of TVD (total-variation-diminishing) Runge- 
Kutta type high order methods introduced in [12]. For example, the third order case is 


( u« = u(°) + AtL(u< 0 >) 
u (2) = Ju<°> + iu^ 1 ) + \AtL{uV) 
U (3) = I u (°) + f u ( 2 ) + |AtL(u( 2 )) 
u (°) = u n , u n+1 = u (3 ) 


(2.9) 
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In the following we will discuss the spatial operator L(u) J in (2.7) and suppress the time 
variable t\ 

(2) The numerical fluxes in (2.7) approximate h ± (x J± i) to high order, where the 
functions h ± (x) are defined implicitly by 

f±(u(x)) = ± jT* h HM < 21 °) 

see [13] for the derivation, r-th order polynomial interpolation is used, based on Newton 
differences. The approximation to h ± (x) is carried out component by component. In the 
following we will use plain letters to represent one component of the corresponding vector in 
bold letters, and we will suppress the superscripts ±. As is pointed out in [13], there is no 
need to construct the function h(x) or its Newton differences explicitly: one can simply use 
the difference tables of /( u(x)). If the undivided differences of /(u(x)) are computed by 

/b')°] = /( u i) (2.11) 

f[j, *] = /[> + 1, * - 1] - /[?, k - 1], k = l, ...,r 
where r is the order of the interpolation polynomial, then the numerical flux f j+ i is obtained 

by 

f, + i = £,C(i-j,k)f[i,k] (2.12) 

3 fc =0 

where i is the left-most point in the stencil used to approximate f j+ 1 , and C(s, k) is defined 
by 

1 a-\-k a-\-k 

°<-*)-(rn)i£ EM (213) 

The small matrix C is independent of the flux, hence can be computed only once and 
then stored. The reader can easily check that (2.12)-(2.13) gives the correct interpolant 
approximation to the function h(x). 

What distinguishes ENO from other finite difference methods is the adaptive stencil 
idea. This idea is realized through the choice of i, the left-most point in the stencil used 
to approximate f j+ i. We start with i = j for computing f+ + i or i = j + 1 for computing 
f~ +i . This is due to upwinding, since the eigenvalues of are all positive, hence the 

information for f + (u) propagates to the right, likewise the information for f (u) propagates 
to the right. The stencil is then expanded point by point, according to the principle of 
choosing the smaller in absolute value of the two relevant differences: 
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if ( abs(f[i , k]).gt.abs(f[i — 1, k ])) i = i — 1 


(2.14) 


for A? = 1, 

A factor can be introduced in (2.14) to bias towards the choice of a linearly stable centered 
stencil [14], This modification can enhance the accuracy in some cases but may cause some 
over-compression effects [14]. Since for nonlinear systems the original ENO scheme works 
well numerically [6], [13], we do not use this modification in this paper. 

The scheme described above is uniformly high order in space and time, measured by local 
truncation error analysis. 

Remark 2.1 Since the schemes described above are conservative, any converged solution 
will be a weak solution of (1.1). It is more difficult to show that the limit solutions are 
admissible. If we take r, the order of the interpolating polynomial, to be zero, and use the 
splitting (2.1 )-(2.2) for a hyperbolic system, we recover the classical Lax- Friedrichs scheme. 
It is well known that the Lax-Friedrichs scheme can be rewritten as a centered scheme 
plus a dissipation term approximately equal to ^aAiu**. If we use the splitting of the 
form (2.3), we can still think about it as a centered scheme plus a dissipation term of the 
form ^aAa:u XI . It is then reasonable, cf. [16], [2], to expect that the scheme converges to 
the admissible weak solutions. Ample numerical tests should be performed to assess the 
convergence and admissability for higher order schemes. The numerical examples in Section 
3 are preliminary results in this direction. 

Remark 2.2 If some fractional step method (e.g. Strang [18]) is used on the splitting 
(2.1), we end up with a scheme of the form (see (2.7)): 

u n+1 = (I + AtL + )(I+AtL-)u n (2.15) 

(for the next time step the two operators may reverse order). For a linear or nonlinear 
problem with smooth solutions, it is very easy to choose stable operators (I + AtL*), due 
to the hyperbolicity of f ± (u) in (2.1). However, this does not necessarily mean that the 
scheme (2.15) is stable, since (I + AiL*) may not commute with each other and may not be 
simultaneously diagonalizable. If the operators satisfy the more restrictive condition 

||I + AiL + || < 1 + 0(A<) ||I + AfL~ || < 1 + 0( Af) (2.16) 

for any consistent norm, the fractional step scheme (2.15) will be stable. 
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3 Numerical Examples 

We use the van der Waals equation (1.2)- (1.3) with RT = 1, a = 0.9 and b = 0.25. The 
graph of the corresponding p(w ) is in Figure 1. The system is elliptic for a < w < f3, 
where a = 0.574912 and /? = 1.036251. The so-called Maxwell line BE in Figure 1, where 
the two shaded areas are equal, intersects the curve of p(w ) at w = m = 0.494273 and 
w = M = 1.405065. The horizontal lines AD and CF in Figure 1 yield 7 = 0.483100 and 

6= 1.918618. 

We use the third order, i.e. (2.9) and (2.12) with r = 3, ENO scheme, described in 
Section 2. If the computational cell is contained completely inside one of the hyperbolic 
regions w < a or w > (3, we use characteristic decompositions (the ENO-LF algorithm 
described in [12], [13]). Otherwise the component by component approximation described 
in Section 2 is used. The splitting used is (2. 3)- (2.5) with a = 2.2. The time step At is 
restricted by a CFL number 0.6, i.e. At < 0.6 M 2 ^) + p{~ t ^ 1 ))^ x where P ( A ) is the 
spectral radius of A. 

All the computations are performed by using a sequence of refined meshes to verify 
convergence, although we typically only show the graphs for one or two fixed meshes. 

We first compute several Riemann problems (1.5): 

(1) {vl,wl) = (l,m), (v R ,w R ) = (1, M) where m and M are the Maxwell values de- 
fined above. This initial condition satisfies the Rankine-Hugoniot condition for a stationary 
jump. Physical principles (Maxwell equal area rule) and many admissibility criteria (e.g. 
the viscosity-capillarity criterion in [15], see (1.4)) indicate that this is an admissible jump. 
Our numerical result shows a stable, sharp jump for this case, Figure 2. We remark that 
here and in what follows, the numerical solution usually has one or two transition points 
in the elliptic region for the phase jump. Apparently they do not cause any trouble to the 
computation. 

(2) {v L ,w L ) = (1,0.54), ( v r ,wr ) = (1,0.54). This initial condition also satisfies the 
Rankine-Hugoniot condition for a stationary jump, but physical principles and many admis- 
sibility criteria (e.g. the viscosity-capillarity criterion in [15], see also [10]) indicate that this 
is not an admissible jump. Our numerical result shows the evolution of this jump into a 
more complicated structure of jumps, Figure 3a, apparently due to the inherent numerical 
viscosity of the scheme (see Remark 2.1). The solution exibits oscillatory behaviors near the 
phase boundary, Figure 3a. This is unpleasant but not surprising since we used component 
by component approximations in cells involving elliptic regions, hence during the process of 
one wave splitting into two or more waves, one or more of them being hyperbolic, oscillations 
occur as a failure of recognition of the corresponding characteristic fields. Similar oscillations 
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also appear for hyperbolic systems if a component by component approximation is used (see, 
e.g. [1]). The oscillations become smaller and more confined when the number of grids is 
increased (Figure 3b), indicating the convergence of the scheme even with these oscillations. 

The background of Figure 3a is computed by the same scheme with 2000 points. It agrees 
with the result with 4000 points hence can be considered as a converged solution. In order 
to check whether this weak solution is admissible as a viscosity-capillarity limit, we plot 
in Figure 3c the numerical solutions of (1.4), for A = i, with e = 0.1,0.01,0.001 and the 
solution of our scheme for (1.2). The solutions to (1.4) are computed by the standard fourth 
order centered scheme with the classical fourth order Runge-Kutta time discretization. We 
verify adequate resolution for the solution of (1.4) for each fixed e by repeatedly refining the 
mesh until the solutions do not change to visual inspection (the largest number of grid points 
used is 8000). Clearly we can see the convergence of the solutions of (1.4) to our solution 
when e — * 0 + in Figure 3c. 

(3) = (1, 0.45), (uH)" 1 ^/?) = (2, 1.5). This case is somewhat easier to compute than 

the previous case, since the initial condition is not a steady nonadmissible weak solution. 
Figure 4a shows the result with 200 grid points on a background of a converged solution 
with 2000 grid points. Figure 4b shows the convergence as e — » 0 + of the solutions of the 
viscosity-capillarity equation (1.4) to our solution. 

We then compute the solutions for smooth initial conditions: 

(1) (v°(x),it; 0 (x)) = (1 — 0.5cos(s), 1 + 0.5.sin(x)). This initial condition crosses the 
elliptic regions (Figure 5). The solution gradually evolves into piecewise smooth solutions 
contained entirely inside one of the two hyperbolic regions w < a and w > /3, connected by 
jumps over the elliptic regions (phase transitions). This seems to agree with the physical 
intuition. 

(2) (u°(x),u; 0 (x)) = (1 — 0.5cos(x), 0.8 + 0.2sin(x)). This initial condition is entirely 
contained in the elliptic region. However, similar to the previous case, the solution gradually 
evolves into piecewise smooth solutions separated by phase transitions, Figure 6. Notice that 
during the evolution oscillations are generated inside the elliptic regions, presumably due to 
the inherent instability of the equation in those regions. These oscillations fade out once the 
solution evolves into the hyperbolic regions. 


4 Concluding Remarks 

Numerical methods for solving systems of conservation laws of mixed hyperbolic-elliptic 
type are investigated, through a flux splitting to write the physical elliptic flux as a sum 
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of two hyperbolic fluxes with positive/negative eigenvalues, then to apply the essentially 
non-oscillatory (ENO) high order finite difference methods on each of them. The method, 
in the simplest first order case, is equivalent to adding a numerical dissipation term with 
a diagonal dissipation matrix. The numerical results on the van der Waals equation of gas 
dynamics indicate that the method can resolve phase boundaries well and can be used as a 
tool to study the evolution of elliptic regions. More numerical tests on different mixed type 
equations constitute current research. 

Acknowledgment: I am grateful to Haitao Fan and Din- Yu Hsieh for bringing my attention 
to the mixed type problems and for many helpful discussions. 
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